Vocal interactions in Thyroptera tricolor

Author

Marcelo Araya-Salas

Published

September 19, 2023

Data analysis for the paper:

Gloriana Chaverri, Maria Sagot, Jennifer L. Stynoski, Marcelo Araya-Salas, Yimen Araya-Ajoy, Martina Nagy, Mirjam Knörnschild, Silvia Chaves-Ramírez, Nicole Rose, Mariela Sánchez-Chavarría, Yelananie Jiménez-Torres, Deilyn Ulloa-Sanabria, Hellen Solís-Hernández, Gerald G. Carter. In review. Contact-calling rates within groups of disc-winged bats vary by caller but not by the receiver’s identity, kinship, or association. Philosophical Transactions of the Royal Society B.

 

Source code and data found at https://github.com/morceglo/Vocal-interactions-Thyroptera-tricolor

Load packages

Code
# load function from sketchy
source("https://raw.githubusercontent.com/maRce10/sketchy/main/R/load_packages.R")

# install/ load packages
sketchy::load_packages(packages = c("igraph", "bipartite", "asnipe",
    "assortnet", "ggplot2", "ggmap", "ecodist", "igraphdata", "statnet",
    "RColorBrewer", "tidyverse", "geosphere", "mapproj", "sna", "stringr",
    "lme4", "glmmTMB", "broom.mixed", "boot", "patchwork", "performance",
    "MASS"))

1 Visualize social networks

1.1 Prepare data

Code
# Using
# https://dshizuka.github.io/networkanalysis/example_make_sparrow_network.html


# Input data
batdat21 = read.csv("./data/raw/associations_2021.csv", header = TRUE,
    sep = ";", colClasses = c("numeric", "character", "Date", "character",
        "character", "character", "character", "character", "character",
        "character", "character", "character", "character"))
batdat22 = read.csv("./data/raw/associations_2022.csv", header = TRUE,
    sep = ";", colClasses = c("numeric", "character", "Date", "character",
        "character", "character", "character", "character", "character",
        "character", "character"))


# Add attribute data
attrib_21 <- read.csv("./data/raw/supp_21.csv", header = TRUE, sep = ";",
    colClasses = c("character", "character"))

attrib_22 <- read.csv("./data/raw/supp_22.csv", header = TRUE, sep = ";",
    colClasses = c("character", "character"))


# Prepare association data 2021
batcols = grep("Bat", colnames(batdat21))
# batcols batdat21[,batcols] gather(batdat21[,batcols])
bat.ids = unique(gather(batdat21[, batcols])$value)
# bat.ids
bat.ids = bat.ids[!is.na(bat.ids)]
# bat.ids
bat.ids_df_21 <- data.frame(bat.ids)
# csv_path <- 'bats_21.csv' write.csv(bat.ids_df_21, csv_path,
# row.names = FALSE)

m1_21 = apply(batdat21[, batcols], 1, function(x) match(bat.ids, x))
# m1_21

m1_21[is.na(m1_21)] = 0
m1_21[m1_21 > 0] = 1
# m1_21

rownames(m1_21) = bat.ids  #rows are bat ids
colnames(m1_21) = paste("group", 1:ncol(m1_21), sep = "_")  #columns are group IDs (just 'group_#')
# m1_21

# rowSums(m1_21) set.seed(2) png(file='./output/Histogram
# captures 2021.png', width=600, height=400) par(oma=c(1,1,1,1))
# # all sides have 3 lines of space par(mar=c(5,4,4,2) + 0.8)
# hist(rowSums(m1_21), main='', xlab='Number of observations',
# ylab='Frequency', cex.lab = 2.5, cex.axis = 1.5) dev.off()

average_n_21 <- mean(rowSums(m1_21))  #average number of times a single bat was observed
# average_n_21

sd_n_21 <- sd(rowSums(m1_21))  #sd number of times a single bat was observed
# sd_n_21

average_gs_21 <- mean(colSums(m1_21))  #average group size
# average_gs_21

sd_gs_21 <- sd(colSums(m1_21))  #sd group size
# sd_gs_21


# Prepare association data 2022
batcols = grep("Bat", colnames(batdat22))
# batcols batdat22[,batcols] gather(batdat22[,batcols])
bat.ids = unique(gather(batdat22[, batcols])$value)
# bat.ids
bat.ids = bat.ids[is.na(bat.ids) == F]
# bat.ids
bat.ids_df_22 <- data.frame(bat.ids)
# csv_path = 'bats_22.csv' write.csv(bat.ids_df_22, csv_path,
# row.names = FALSE)

m1_22 = apply(batdat22[, batcols], 1, function(x) match(bat.ids, x))
# m1_22

m1_22[is.na(m1_22)] = 0
m1_22[m1_22 > 0] = 1
# m1_22

rownames(m1_22) = bat.ids  #rows are bat ids
colnames(m1_22) = paste("group", 1:ncol(m1_22), sep = "_")  #columns are group IDs (just 'group_#')
# m1_22

# rowSums(m1_22) set.seed(2) png(file='./output/Histogram
# captures 2022.png', width=600, height=400) par(oma=c(1,1,1,1))
# # all sides have 3 lines of space par(mar=c(5,4,4,2) + 0.8)
# hist (rowSums(m1_22), main='', xlab='Number of observations',
# ylab='Frequency', cex.lab = 2.5, cex.axis = 1.5) dev.off()

average_n_22 <- mean(rowSums(m1_22))  #average number of times a single bat was observed
# average_n_22

sd_n_22 <- sd(rowSums(m1_22))  #sd number of times a single bat was observed
# sd_n_22

average_gs_22 <- mean(colSums(m1_22))  #average group size
# average_gs_22

sd_gs_22 <- sd(colSums(m1_22))  #average group size
# sd_gs_22

1.2 Plot networks

Code
#Create adjacency data for associations (2021)
adj_21=get_network(t(m1_21), data_format="GBI", association_index = "SRI") # the adjacency matrix
Generating  76  x  76  matrix
Code
g_21=graph_from_adjacency_matrix(adj_21, "undirected", weighted=T) #the igraph object
V(g_21)$sex <- factor(attrib_21[match(V(g_21)$name, attrib_21$bat_id), "sex"])

mismatched_names <- setdiff(V(g_21)$name, attrib_21$bat_id) #check for issues


#Create adjacency data for associations (2022)
adj_22=get_network(t(m1_22), data_format="GBI", association_index = "SRI") # the adjacency matrix
Generating  76  x  76  matrix
Code
g_22=graph_from_adjacency_matrix(adj_22, "undirected", weighted=T) #the igraph object
V(g_22)$sex <- factor(attrib_22[match(V(g_22)$name, attrib_22$bat_id), "sex"])

mismatched_names <- setdiff(V(g_22)$name, attrib_22$bat_id) #check for issues


#Select individuals for which more than 10 observations were made (not incorporated into adjacency data for associations, above)
#m2_21=m1_21[which(rowSums(m1_21)>10),]
#adj_21=get_network(t(m2_21), data_format="GBI","SRI")
#write.csv(adj_21, "assoc_21.csv", row.names=TRUE)
#rowSums(m2_21)

#m2_22=m1_22[which(rowSums(m1_22)>10),]
#adj_22=get_network(t(m2_22), data_format="GBI","SRI")
#rowSums(m2_22)
#write.csv(adj_22, "assoc_22.csv", row.names=TRUE)
#rowSums(m2_22)


#Create communities (2021)
com_21=fastgreedy.community(g_21) #community detection method
# com_21[[11]]

# Assign community names
mapping <- c(
  11,  # Old membership 1 should be corrected to 11
  8,   # Old membership 2 should be corrected to 8
  1,   # Old membership 3 should be corrected to 1
  2,   # Old membership 4 should be corrected to 2
  6,   # Old membership 5 should be corrected to 6
  9,   # Old membership 6 should be corrected to 9
  4,   # Old membership 7 should be corrected to 4
  10,  # Old membership 8 should be corrected to 10
  7,   # Old membership 9 should be corrected to 7
  5,   # Old membership 10 should be corrected to 5
  3    # Old membership 11 should be corrected to 3
)

corrected_memberships <- mapping[com_21$membership]

com_21$membership <- corrected_memberships

community_assignments <- com_21$membership

color_mapping <- c(
  "female" = "slateblue",
  "male" = "gold",
  "unknown" = "gray"
)
node_sex <- V(g_21)$sex
node_colors <- sapply(node_sex, function(sex) {
  color_mapping[sex]
})
V(g_21)$color <- node_colors

#Prepare to save figure 
# png(file="./output/Association networks 2021.png",
    # width=600, height=600)
 par(oma=c(1,1,1,1), mar=c(5,4, 4, 2) + 0.8, mfrow = c(1, 2))
# Iterate through com_21 and assign community names to nodes in g_21
plot(
  g_21, 
  vertex.color = node_colors,  # Use the calculated node colors
  vertex.size = 10,  # Adjust the size of the nodes as needed
  vertex.label = community_assignments, # Use the extracted community assignments as labels
  edge.width=E(g_21)$weight*10
)

title(main = "2021", cex.main = 2)

legend("topleft", legend = names(color_mapping), fill = color_mapping, cex = 1.5)

#Create communities (2022) 
com_22=fastgreedy.community(g_22) #community detection method

# Assign community names
mapping <- c(
  3,  
  4, 
  1,  
  9, 
  10, 
  11,  
  7, 
  5, 
  6, 
  2, 
  12,
  8
)

corrected_memberships <- mapping[com_22$membership]

com_22$membership <- corrected_memberships

community_assignments <- com_22$membership

color_mapping <- c(
  "female" = "slateblue",
  "male" = "gold",
  "unknown" = "gray"
)
node_sex <- V(g_22)$sex
node_colors <- sapply(node_sex, function(sex) {
  color_mapping[sex]
})
V(g_22)$color <- node_colors

#Prepare to save figure 
# png(file="./output/Association networks 2022.png",width=600, height=600)

# Iterate through com_21 and assign community names to nodes in g_21
plot(
  g_22, 
  vertex.color = node_colors,  # Use the calculated node colors
  vertex.size = 10,  # Adjust the size of the nodes as needed
  vertex.label = community_assignments, # Use the extracted community assignments as labels
  edge.width=E(g_22)$weight*10
)

title(main = "2022", cex.main = 2)

1.3 Estimate modulatity

Modularity:

  • 2021: 0.87
  • 2022: 0.9

2 Calling rate repeatability

2.1 Negative binomial model

Code
# Data manipulation Joining the two data sets
d1 <- read.csv("./data/raw/vocal_interactions_2021.csv", sep = ";")
d2 <- read.csv("./data/raw/vocal_interactions_2022.csv", sep = ";")
colnames(d2)[6] <- "Dyad"
d <- rbind(d1, d2)

## Adding observation level random effect
d$ob <- 1:nrow(d)

# Stats models Model for inquiry
mod_inquiry <- glmer.nb(n_inquiry ~ 1 + (1 | Bat_roosting) + (1 |
    Bat_flying) + (1 | Group), data = d)
summary(mod_inquiry)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(2.7219)  ( log )
Formula: n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |  
    Group)
   Data: d

     AIC      BIC   logLik deviance df.resid 
  4596.7   4617.2  -2293.3   4586.7      446 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.62529 -0.48778  0.07217  0.47606  2.82516 

Random effects:
 Groups       Name        Variance  Std.Dev. 
 Bat_roosting (Intercept) 8.029e-11 8.961e-06
 Bat_flying   (Intercept) 2.348e-01 4.845e-01
 Group        (Intercept) 2.141e-02 1.463e-01
Number of obs: 451, groups:  Bat_roosting, 74; Bat_flying, 73; Group, 14

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.12062    0.07845   52.53   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Code
### Extracting variance components
mu_l <- fixef(mod_inquiry)[1]
mu <- exp(mu_l)
theta <- getME(mod_inquiry, "glmer.nb.theta")

VBF_I_l <- VarCorr(mod_inquiry)$Bat_flying[1, 1]  #Variance Bats Flying
VBF_I <- (exp(VBF_I_l) - 1) * exp(2 * mu_l + VBF_I_l)

VBR_I_l <- VarCorr(mod_inquiry)$Bat_roosting[1, 1]  #Variance Bats Roosting
VBR_I <- (exp(VBR_I_l) - 1) * exp(2 * mu_l + VBR_I_l)

VBG_I_l <- VarCorr(mod_inquiry)$Group[1, 1]  #Variance GRoup
VBG_I <- (exp(VBG_I_l) - 1) * exp(2 * mu_l + VBG_I_l)


NBV_I <- mu + mu^2/theta

#### Repeatabilty roosting individual
VBR_I/(VBR_I + VBF_I + NBV_I + VBG_I)
 (Intercept) 
1.084486e-10 
Code
#### Repeatabilty flying individual
VBF_I/(VBR_I + VBF_I + NBV_I)
(Intercept) 
   0.465914 
Code
# Stats models Model for inquiry2
d2 <- d[d$n_response > 0, ]
mod_inquiry <- glmer.nb(n_inquiry ~ 1 + (1 | Bat_roosting) + (1 |
    Bat_flying) + (1 | Group), data = d2)
summary(mod_inquiry)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(4.214)  ( log )
Formula: n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |  
    Group)
   Data: d2

     AIC      BIC   logLik deviance df.resid 
  2606.7   2624.4  -1298.4   2596.7      250 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.93767 -0.52024  0.04805  0.49889  2.73912 

Random effects:
 Groups       Name        Variance  Std.Dev. 
 Bat_roosting (Intercept) 7.930e-13 8.905e-07
 Bat_flying   (Intercept) 2.221e-01 4.712e-01
 Group        (Intercept) 3.879e-02 1.970e-01
Number of obs: 255, groups:  Bat_roosting, 71; Bat_flying, 69; Group, 14

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.21595    0.08999   46.85   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Code
### Extracting variance components
mu_l <- fixef(mod_inquiry)[1]
mu <- exp(mu_l)
theta <- getME(mod_inquiry, "glmer.nb.theta")

VBF_I_l <- VarCorr(mod_inquiry)$Bat_flying[1, 1]  #Variance Bats Flying
VBF_I <- (exp(VBF_I_l) - 1) * exp(2 * mu_l + VBF_I_l)

VBR_I_l <- VarCorr(mod_inquiry)$Bat_roosting[1, 1]  #Variance Bats Roosting
VBR_I <- (exp(VBR_I_l) - 1) * exp(2 * mu_l + VBR_I_l)

VBG_I_l <- VarCorr(mod_inquiry)$Group[1, 1]  #Variance GRoup
VBG_I <- (exp(VBG_I_l) - 1) * exp(2 * mu_l + VBG_I_l)


NBV_I <- mu + mu^2/theta

#### Repeatabilty roosting individual
VBR_I/(VBR_I + VBF_I + NBV_I)
 (Intercept) 
1.409571e-12 
Code
#### Repeatabilty flying individual
VBF_I/(VBR_I + VBF_I + NBV_I)
(Intercept) 
  0.5519081 
Code
## Model for response
mod_response <- glmer.nb(n_response ~ 1 + (1 | Bat_roosting) + (1 |
    Bat_flying) + (1 | Group), data = d)
summary(mod_response)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(0.1622)  ( log )
Formula: n_response ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |  
    Group)
   Data: d

     AIC      BIC   logLik deviance df.resid 
  3742.3   3762.8  -1866.1   3732.3      446 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.4026 -0.4017 -0.3616  0.2471  4.1603 

Random effects:
 Groups       Name        Variance  Std.Dev. 
 Bat_roosting (Intercept) 1.118e+00 1.057e+00
 Bat_flying   (Intercept) 5.331e-13 7.301e-07
 Group        (Intercept) 1.489e-12 1.220e-06
Number of obs: 451, groups:  Bat_roosting, 74; Bat_flying, 73; Group, 14

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)    4.074      0.213   19.12   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Code
### Extracting variance components
mu_l <- fixef(mod_response)[1]  #+ fixef(mod_response)[2]*mean(d$n_inquiry, na.rm=TRUE)
mu <- exp(mu_l)
theta <- getME(mod_response, "glmer.nb.theta")

VBF_R_l <- VarCorr(mod_response)$Bat_flying[1, 1]  #Variance Bats Flying
VBF_R <- (exp(VBF_R_l) - 1) * exp(2 * mu_l + VBF_R_l)

VBR_R_l <- VarCorr(mod_response)$Bat_roosting[1, 1]  #Variance Bats Roosting
VBR_R <- (exp(VBR_R_l) - 1) * exp(2 * mu_l + VBR_R_l)

VBG_R_l <- VarCorr(mod_response)$Group[1, 1]  #Variance GRoup
VBG_R <- (exp(VBG_R_l) - 1) * exp(2 * mu_l + VBG_R_l)


NBV_R <- mu + mu^2/theta

#### Repeatabilty roosting individual
VBR_R/(VBR_R + VBF_R + NBV_R + VBG_R)
(Intercept) 
     0.5044 
Code
#### Repeatabilty flying individual
VBF_R/(VBR_R + VBF_R + NBV_R + VBG_R)
 (Intercept) 
4.274015e-14 

2.2 Poisson model

Code
# #Data manipulation ## Joining the two data sets
# d1<-read.csv('vocal_interactions_2021.csv', sep=';')
# d2<-read.csv('vocal_interactions_2022.csv', sep=';')
# colnames(d2)[6]<-'Dyad' d<-rbind(d1,d2)

## Adding observation level random effect
d$ob <- 1:nrow(d)

# Stats models Model for inquiry with all observations
mod_inquiry <- glmer(n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) +
    (1 | Group) + (1 | ob), data = d, family = "poisson")

### Extracting variance components
VBF_I <- VarCorr(mod_inquiry)$Bat_flying[1, 1]  #Variance Bats Flying
VBR_I <- VarCorr(mod_inquiry)$Bat_roosting[1, 1]  #Variance Bats Roosting
VBG_I <- VarCorr(mod_inquiry)$Group[1, 1]  #Variance Group
VOD_I <- VarCorr(mod_inquiry)$ob[1, 1]  #Variance Overddispersion

PV_I <- log(1/exp(fixef(mod_inquiry)[1]) + 1)  #Variance poisson processs

#### Repeatabilty roosting individual
VBR_I/(VBR_I + VOD_I + VBF_I + PV_I + VBG_I)
(Intercept) 
 0.07582562 
Code
#### Repeatabilty flying individual
VBF_I/(VBR_I + VOD_I + VBF_I + PV_I + VBG_I)
(Intercept) 
  0.3913205 
Code
## Model for inquiry only when someone responded
d2 <- d[d$n_response > 0, ]
mod_inquiry2 <- glmer(n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) +
    (1 | ob), data = d2, family = "poisson")

### Extracting variance components
VBF_I <- VarCorr(mod_inquiry2)$Bat_flying[1, 1]  #Variance Bats Flying
VBR_I <- VarCorr(mod_inquiry2)$Bat_roosting[1, 1]  #Variance Bats Roosting
VBG_I <- VarCorr(mod_inquiry)$Group[1, 1]
VOD_I <- VarCorr(mod_inquiry2)$ob[1, 1]  #Variance Overddispersion
PV_I <- log(1/exp(fixef(mod_inquiry2)[1]) + 1)  #Variance poisson processs

#### Repeatabilty roosting individual
VBR_I/(VBR_I + VOD_I + VBF_I + PV_I + VBG_I)
(Intercept) 
 0.02516122 
Code
#### Repeatabilty flying individual
VBF_I/(VBR_I + VOD_I + VBF_I + PV_I + VBG_I)
(Intercept) 
  0.5240666 
Code
## Model for response
mod_response <- glmer(n_response ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) +
    (1 | Group) + (1 | ob), data = d, family = "poisson")
summary(mod_response)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: n_response ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |  
    Group) + (1 | ob)
   Data: d

     AIC      BIC   logLik deviance df.resid 
  3816.4   3837.0  -1903.2   3806.4      446 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-0.75743 -0.39416 -0.00652  0.01611  0.17248 

Random effects:
 Groups       Name        Variance Std.Dev.
 ob           (Intercept) 8.3383   2.8876  
 Bat_roosting (Intercept) 7.8244   2.7972  
 Bat_flying   (Intercept) 0.4551   0.6746  
 Group        (Intercept) 1.0819   1.0401  
Number of obs: 451, groups:  
ob, 451; Bat_roosting, 74; Bat_flying, 73; Group, 14

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 1.1595061  0.0004604    2518   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.15353 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?
Code
### Extracting variance components
VBF_R <- VarCorr(mod_response)$Bat_flying[1, 1]  #Variance Bats Flying
VBR_R <- VarCorr(mod_response)$Bat_roosting[1, 1]  #Variance Bats Roosting
VBG_R <- VarCorr(mod_response)$Group[1, 1]
VOD_R <- VarCorr(mod_response)$ob[1, 1]  #Variance Overddispersion
PV_R <- log(1/exp(fixef(mod_response)[1]) + 1)  #Variance poisson processs

## Repeatabilty roosting individual
VBR_R/(VBR_R + VOD_R + VBF_R + PV_R + VBG_R)
(Intercept) 
  0.4353549 
Code
## Repeatabilty flying individual
VBF_R/(VBR_R + VOD_R + VBF_R + PV_R + VBG_R)
(Intercept) 
 0.02532077 

3 Calling, association and kinship

3.1 Prepare data

Code
# Code by Gerald Carter, gcarter1640@gmail.com

### clean and wrangle data

# create function to convert matrix to data-frame
matrix_to_df <- function(m1){
  data.frame(dyad = paste(rownames(m1)[col(m1)], colnames(m1)[row(m1)], sep="_"),
             value = c(t(m1)), stringsAsFactors = FALSE)
}

# get kinship data
k <- as.matrix(read.csv('./data/raw/KinshipMatrix.csv', sep= ",", row.names = 1))

# check symmetry
mean(k==t(k), na.rm=T)

# get sex data
sex1 <- read.csv('supp_21.csv', sep= ";") 
sex2 <- read.csv('supp_22.csv', sep= ";") 

# get association data
a21 <- read.csv('associations_2021.csv', sep= ";")
a22 <- read.csv('associations_2022.csv', sep= ";")

# get calling data
i21 <- read.csv('vocal_interactions_2021.csv', sep= ";")
i22 <- read.csv('vocal_interactions_2022.csv', sep= ";")

# tidy sex data
sex <- 
  rbind(sex1,sex2) %>% 
  mutate(bat= paste0('bat', bat_id)) %>% 
  group_by(bat, sex) %>% 
  summarize(n= n()) %>% 
  dplyr::select(bat, sex) %>% 
  ungroup()

# tidy 2021 association data
assoc21 <- 
  a21 %>% 
  pivot_longer(Bat1:Bat10, names_to = 'names', values_to = "bat") %>% 
  filter(!is.na(bat)) %>% 
  mutate(bat= paste0('bat', bat)) %>% 
  mutate(group= paste0('group', Group)) %>%  
  mutate(date= dmy(Date)) %>% 
  dplyr::select(obs= data_id, group, date, bat)

# tidy 2022 association data
assoc22 <- 
  a22 %>% 
  pivot_longer(Bat1:Bat8, names_to = 'names', values_to = "bat") %>% 
  filter(!is.na(bat)) %>% 
  mutate(bat= paste0('bat', bat)) %>% 
  mutate(group= paste0('group', Group)) %>%  
  mutate(date= dmy(Date)) %>% 
  dplyr::select(obs= data_id, group, date, bat)

# combine association data
assoc <- 
  rbind(assoc21, assoc22) %>% 
  mutate(obs= paste(obs,group,date, sep= "_")) %>% 
  dplyr::select(bat, obs) %>% 
  as.data.frame()

# get number of observations per bat
obs_per_bat <- 
  assoc %>% 
  group_by(bat) %>% 
  summarize(n.obs=n()) %>% 
  arrange(desc(n.obs))

# range is 1 to 51 times
range(obs_per_bat$n.obs)

# pick threshold of observations for including bats in analysis
# if they have fewer observations than we make their association rates NA
threshold <- 4

# plot number of observations per bat
assoc %>% 
  group_by(bat) %>% 
  summarize(n.obs=n()) %>% 
  ggplot(aes(x=n.obs))+
    geom_histogram(fill="light blue", color="black")+
    geom_vline(xintercept= threshold, color= 'red', size=1)+
    xlab("number of observations of bat")+
    ylab("count of bats")

# get bats seen fewer than threshold number
low.n.bats <- 
  assoc %>% 
  group_by(bat) %>% 
  summarize(n.obs=n()) %>% 
  filter(n.obs < threshold) %>% 
  pull(bat)
low.n.bats
length(low.n.bats)
# 14 bats were seen fewer than 4 times

### get SRI from asnipe

# get association rates as group-by-individual
gbi <- get_group_by_individual(assoc, data_format="individuals")

# get association network of SRI
assoc.net<- get_network(association_data=gbi, data_format = "GBI", association_index = "SRI")

# check symmetry
mean(assoc.net==t(assoc.net))

# get SRI values for every undirected pair as dataframe
SRI <- matrix_to_df(assoc.net) 

# tidy kinship
kinship <- 
  k %>% 
  matrix_to_df() %>% 
  separate(dyad, into= c("bat1", "bat2")) %>% 
  mutate(bat2= str_remove(bat2, "X")) %>% 
  mutate(bat1= paste0('bat', bat1), bat2= paste0('bat', bat2)) %>% 
  filter(bat1!=bat2) %>% 
  # label undirected pair
  mutate(dyad= ifelse(bat1<bat2, paste(bat1,bat2, sep="_"), paste(bat2,bat1, sep="_"))) %>% 
  group_by(dyad) %>% 
  summarize(kinship= first(value)) 

# count trials
nrow(i21)
nrow(i22)

# trials with missing data
rbind(i21,i22) %>% 
  as_tibble() %>% 
  filter(is.na(n_response)) %>% 
  nrow()
#2

# trials where neither bat called
rbind(i21,i22) %>% 
  as_tibble() %>% 
  filter(n_inquiry==0) %>% 
  nrow()
# 5

# fix and tidy calling data
calling <- 
  rbind(i21,i22) %>% 
  mutate(Date= as.character(mdy(Date))) %>% 
  mutate(group= paste(substr(Date, start=1, stop=4),Group, sep="_")) %>% 
  # need to recreate these labels because messed up by excel in raw data
  separate(Dyad, into=c("bat_flying", "bat_roosting"), remove=F) %>% 
  mutate(bat_flying= paste0('bat', bat_flying), bat_roosting= paste0('bat', bat_roosting)) %>% 
  # label undirected dyads
  mutate(dyad= ifelse(bat_flying<bat_roosting, paste(bat_flying, bat_roosting, sep= "_"), paste(bat_roosting, bat_flying, sep= "_"))) %>% 
  dplyr::select(date= Date, group, bat_flying, bat_roosting, dyad, flight_time, n_inquiry, n_response)

# combine all data
d <- 
  # start with calling data
  calling %>% 
  # add SRI 
  mutate(sri= SRI$value[match(.$dyad, SRI$dyad)]) %>% 
  # add kinship
  mutate(kinship= kinship$kinship[match(.$dyad, kinship$dyad)]) %>% 
  # add sexes of flying and roosting bats
  mutate(flying.sex= sex$sex[match(.$bat_flying, sex$bat)]) %>% 
  mutate(roosting.sex= sex$sex[match(.$bat_roosting, sex$bat)]) %>% 
  # label sex combination for flying-->roosting dyad
  mutate(dyad.sex= paste(flying.sex, roosting.sex, sep= "-->")) %>% 
  # label sex combination for dyad (male, female, both)
  mutate(udyad.sex = case_when(
    flying.sex == 'female' & roosting.sex == 'female' ~ "female-female",
    flying.sex == 'male' & roosting.sex == 'male' ~ "male-male",
    TRUE ~ "mixed-sex")) %>% 
  # replace zero sri (never previously seen together) with NA because association is not 0
  mutate(sri = if_else(sri==0, NA, sri)) %>% 
  # if flying bat seen fewer than 4 times, then SRI is NA
  mutate(sri = if_else(bat_flying %in% low.n.bats, NA, sri)) %>% 
  # if roosting bat seen fewer than 4 times, then SRI is NA
  mutate(sri = if_else(bat_roosting %in% low.n.bats, NA, sri))
           
  
# inspect and fix cases where flying bats in more than 2 groups
t <- 
  d %>% 
  group_by(bat_flying, date, group) %>%
  summarize(n=n()) %>% 
  arrange(date) %>% 
  group_by(bat_flying, group) %>% 
  summarize(n=n()) %>%   
  group_by(bat_flying) %>% 
  summarize(n=n()) %>% 
  filter(n>2) %>% 
  pull(bat_flying)
d %>% 
  dplyr::select(date, group, bat_flying) %>% 
  filter(bat_flying %in% t) %>% 
  arrange(bat_flying) %>% 
  as.data.frame()
# some of these seem impossible and must be errors

# fix impossible group labels
d$group[which(d$date=="2021-07-03" & d$bat_flying=="bat982126051278521")] <- "2021_9"
d$group[which(d$date=="2021-07-03" & d$bat_flying=="bat982126058484300")] <- "2021_9"

# fix cases where roosting bats in more than 2 groups
t <- 
  d %>% 
  group_by(bat_roosting, date, group) %>%
  summarize(n=n()) %>% 
  arrange(date) %>% 
  group_by(bat_roosting, group) %>% 
  summarize(n=n()) %>%   
  group_by(bat_roosting) %>% 
  summarize(n=n()) %>% 
  filter(n>2) %>% 
  pull(bat_roosting)
d %>% 
  dplyr::select(date, group, bat_roosting) %>% 
  filter(bat_roosting %in% t) %>% 
  arrange(bat_roosting) %>% 
  as.data.frame()

# fix group labels
d$group[which(d$date=="2021-07-03" & d$bat_roosting=="bat982126058484300")] <- "2021_9"

# group 9 split into groups 9,1 and 9,2 during 2022
# There are interesting movements between group 9 and 10 and between groups 9,1 and 9,2
# For this analysis, I'm going to consider group 9,1 and 9,2 as the same group
d$group[which(d$group=="2022_9,1" | d$group=="2022_9,2")] <- "2021_9"  

# remove trials without inquiry or response calls
d <- 
  d %>% # 453 rows
  as_tibble() %>% 
  filter(! (is.na(n_inquiry) & is.na(n_response))) %>%   # 451 rows
  filter(n_inquiry>0) %>% # 446 rows
  mutate(year= substr(date, 1,4)) # add year

# save data as csv
write.csv(d, file = paste0('./processed/calling-association-kinship-data.csv'), row.names= F)

3.2 Explore data

Code
# set ggplot theme
theme_set(theme_bw(base_size = 14))

# get data
raw <- read.csv("./data/processed/calling-association-kinship-data.csv")


###### Create functions to estimate confindence intervals
###### function to get mean and 95% CI of values x via
###### bootstrapping
boot_ci <- function(x, perms = 5000, bca = F) {
    get_mean <- function(x, d) {
        return(mean(x[d]))
    }
    x <- as.vector(na.omit(x))
    mean <- mean(x)
    if (bca) {
        boot <- boot.ci(boot(data = x, statistic = get_mean, R = perms,
            parallel = "multicore", ncpus = 4), type = "bca")
        low <- boot$bca[1, 4]
        high <- boot$bca[1, 5]
    } else {
        boot <- boot.ci(boot(data = x, statistic = get_mean, R = perms,
            parallel = "multicore", ncpus = 4), type = "perc")
        low <- boot$perc[1, 4]
        high <- boot$perc[1, 5]
    }
    c(low = low, mean = mean, high = high, N = round(length(x)))
}


# function to get mean and 95% CI via bootstrapping of values y
# within grouping variable x
boot_ci2 <- function(d = d, y = d$y, x = d$x, perms = 5000, bca = F) {
    df <- data.frame(effect = unique(x))
    df$low <- NA
    df$mean <- NA
    df$high <- NA
    df$n.obs <- NA
    for (i in 1:nrow(df)) {
        ys <- y[which(x == df$effect[i])]
        if (length(ys) > 1 & var(ys) > 0) {
            b <- boot_ci(y[which(x == df$effect[i])], perms = perms,
                bca = bca)
            df$low[i] <- b[1]
            df$mean[i] <- b[2]
            df$high[i] <- b[3]
            df$n.obs[i] <- b[4]
        } else {
            df$low[i] <- min(ys)
            df$mean[i] <- mean(ys)
            df$high[i] <- max(ys)
            df$n.obs[i] <- length(ys)
        }
    }
    df
}

# look at data head(raw) group is the year and social group dyad
# is the undirected pair (flying and roosting bat) in
# alphanumeric order sri is simple ratio index of association
# flying.sex is sex of flying bat roosting.sex is sex of
# roosting bat dyad.sex is the sex of the flying and roosting
# bat udyad.sex is females, males, or mixed

# add a few variables to the data
d <- raw %>%
    # get log count of inquiry calls
mutate(log_inquiry = log(n_inquiry)) %>%
    # rescale kinship so 1 unit is 0.5
mutate(kinship2 = kinship * 2) %>%
    # if the roosting bat leaves the roost (escapes) then flight
    # time is < 300 s
mutate(escape = as.numeric(flight_time < 300)) %>%
    # convert flight time from seconds to minutes (for easier
    # interpretation)
mutate(flight_time2 = flight_time/60)


###### Describe the data-------------

# how many trials?  d %>% nrow() 446

# # how many undirected pairs have vocal data?  d %>% pull(dyad)
# %>% unique() %>% length 139 undirected pairs

# d %>% group_by(bat_flying, bat_roosting) %>% summarize(n=n())
# 254 directed pairs

# how many groups?  n_distinct(d$group) 23

# how many of these have association data d %>% group_by(dyad)
# %>% summarize(sri= mean(sri)) %>% filter(!is.na(sri)) 133
# pairs have association

# how many of these have kinship data d %>% group_by(dyad) %>%
# summarize(kinship= mean(kinship)) %>% filter(!is.na(kinship))
# 82 pairs have kinship data

# how many related undirected pairs d %>% group_by(dyad) %>%
# summarize(kinship= mean(kinship)) %>% filter(kinship>0) 32 kin
# pairs

# how many unrelated undirected pairs d %>% group_by(dyad) %>%
# summarize(kinship= mean(kinship)) %>% filter(kinship==0) 50
# nonkin pairs

# what is mean kinship in group?
grp.kin <- d %>%
    group_by(dyad, group) %>%
    summarize(kinship = mean(kinship, na.rm = T)) %>%
    group_by(group) %>%
    summarize(kinship = mean(kinship, na.rm = T), n = n()) %>%
    as.data.frame()

set.seed(123)
# grp.kin %>% pull(kinship) %>% boot_ci(bca=T) 0.23, 95% CI =
# [0.16, 0.31]

# how many groups have kin?  grp.kin %>% filter(kinship>0) 17

# how many groups have no kin?  grp.kin %>% filter(kinship==0) 4

# plot distribution of kinship
d %>%
    ggplot(aes(x = kinship)) + facet_wrap(~udyad.sex, ncol = 1, scales = "free_y") +
    geom_histogram(fill = "light blue", color = "black") + ggtitle("pairwise kinship")

Code
# # check categories d %>% filter(!is.na(kinship)) %>%
# group_by(dyad) %>% summarize(kinship= mean(kinship)) %>%
# group_by(kinship) %>% summarize(n=n())

##### Effect of kinship on association-------------

# plot distribution of assoc
d %>%
    ggplot(aes(x = sri)) + facet_wrap(~udyad.sex, ncol = 1, scale = "free_y") +
    geom_histogram(fill = "light blue", color = "black") + xlab("association rate") +
    ggtitle("pairwise SRI values")

Code
# looks ok

# what is mean and 95% CI of assoc among flying and roosting
# bats?  set.seed(123) d %>% group_by(dyad) %>% summarize(sri =
# mean(sri)) %>% pull(sri) %>% boot_ci(bca=T) 0.51, 95% CI =
# [0.47, 0.55] a bit low, expected from past work is 0.76

# now only nonkin
set.seed(123)
k1 <- d %>%
    filter(kinship == 0) %>%
    group_by(dyad) %>%
    summarize(sri = mean(sri)) %>%
    pull(sri) %>%
    boot_ci(bca = T) %>%
    c(kinship = 0)
# k1 0.572, 95% CI = [0.500, 0.636]

# now only kin set.seed(123) d %>% filter(kinship>0) %>%
# group_by(dyad) %>% summarize(sri = mean(sri)) %>% pull(sri)
# %>% boot_ci(bca=T) 0.686, 95% CI = [0.631, 0.741]

# now only close kin
set.seed(123)
k2 <- d %>%
    filter(kinship == 0.5) %>%
    group_by(dyad) %>%
    summarize(sri = mean(sri)) %>%
    pull(sri) %>%
    boot_ci(bca = T) %>%
    c(kinship = 0.5)
# k2 0.687, 95% CI = [0.628, 0.750]

# now only 0.25 kin
set.seed(123)
k3 <- d %>%
    filter(kinship == 0.25) %>%
    group_by(dyad) %>%
    summarize(sri = mean(sri)) %>%
    pull(sri) %>%
    boot_ci(bca = T) %>%
    c(kinship = 0.25)
# k3 0.684, 95% CI = [0.542, 0.783]

# compile means
k <- rbind(k1, k2, k3) %>%
    data.frame() %>%
    mutate(kin = as.character(kinship)) %>%
    mutate(kin2 = kinship > 0) %>%
    mutate(assoc = mean)

# plot association by kinship
(kin.dyads.plot <- d %>%
    filter(!is.na(kinship)) %>%
    group_by(dyad) %>%
    summarize(kinship = mean(kinship), assoc = mean(sri), udyad.sex = first(udyad.sex)) %>%
    mutate(kin2 = kinship > 0) %>%
    mutate(kin = as.character(kinship)) %>%
    ggplot(aes(x = kin, y = assoc, group = kin, color = kin2)) + geom_jitter(size = 2,
    width = 0.1, height = 0, aes(shape = udyad.sex)) + geom_point(data = k,
    position = position_nudge(x = 0.25), size = 3) + geom_errorbar(data = k,
    aes(ymin = low, ymax = high, width = 0.1), position = position_nudge(x = 0.25),
    size = 1) + xlab("kinship") + ylab("association rate (simple ratio index)") +
    scale_color_manual(values = c("darkgrey", "darkred")) + theme(legend.position = "none",
    legend.title = element_blank()))

Code
# save as PDF ggsave( 'kin_association.pdf', plot =
# kin.dyads.plot, scale = 1, width = 3, height = 5, units =
# 'in', dpi = 300)


##### Effect of flight time on calling ----------

### How does flight time and count of inquiry calls predict
### count of response calls?

# plot number of response calls by flight time
t1 <- d %>%
    mutate(kinship = kinship > 0.1) %>%
    ggplot(aes(x = flight_time, y = n_response)) + geom_point(size = 2,
    alpha = 0.3, aes(color = kinship, shape = kinship)) + geom_smooth(method = "glm.nb") +
    scale_color_manual(values = c("darkgrey", "darkred")) + xlab("flight time (seconds)") +
    ylab("response call count") + theme(legend.position = "none")

# plot number of inquiry calls by flight time
t2 <- d %>%
    mutate(kinship = kinship > 0.1) %>%
    ggplot(aes(x = flight_time, y = n_inquiry)) + geom_point(size = 2,
    alpha = 0.3, aes(color = kinship, shape = kinship)) + geom_smooth(method = "glm.nb") +
    scale_color_manual(values = c("darkgrey", "darkred")) + xlab("flight time (seconds)") +
    ylab("inquiry call count") + theme(legend.position = "none")

# plot together
t <- t1 + t2
# ggsave( filename= 'flight_time.pdf', plot = t, scale = 1,
# width = 7, height = 7, units = c('in', 'cm', 'mm', 'px'), dpi
# = 300)

3.3 Statistical analysis

Code
# try poisson model for effect of flight time on inquiry calling
t <- glmer(n_inquiry ~ flight_time2 + (1 | bat_flying) + (1 | bat_roosting) +
    (1 | group), data = d, family = poisson)
check_overdispersion(t)
# Overdispersion test

       dispersion ratio =    6.288
  Pearson's Chi-Squared = 2773.032
                p-value =  < 0.001
Code
# negative binomial model (NBM)
t <- glmmTMB(n_inquiry ~ flight_time2 + (1 | bat_flying) + (1 | bat_roosting) +
    (1 | group), data = d, ziformula = ~0, family = nbinom2)
check_overdispersion(t)
# Overdispersion test

       dispersion ratio =   0.786
  Pearson's Chi-Squared = 345.749
                p-value =       1
Code
t2 <- tidy(t, conf.int = TRUE, exponentiate = T, effects = "fixed",
    conf.method = "profile")
t2
# A tibble: 2 × 9
  effect component term         estimate std.error statistic  p.value conf.low
  <chr>  <chr>     <chr>           <dbl>     <dbl>     <dbl>    <dbl>    <dbl>
1 fixed  cond      (Intercept)     11.9     1.51        19.5 7.77e-85     9.30
2 fixed  cond      flight_time2     1.43    0.0336      15.1 2.60e-51     1.36
# ℹ 1 more variable: conf.high <dbl>
Code
# for every minute of flight time, the inquiry call count
# increases by a factor of 1.43 [1.36, 1.49]

# try poisson model for effect of flight time on response
# calling
t <- glmer(n_response ~ flight_time2 + log_inquiry + (1 | bat_flying) +
    (1 | bat_roosting) + (1 | group), data = d, family = poisson)
check_overdispersion(t)
# Overdispersion test

       dispersion ratio =    63.477
  Pearson's Chi-Squared = 27930.076
                p-value =   < 0.001
Code
# try NBM
t <- glmmTMB(n_response ~ flight_time2 + log_inquiry + (1 | bat_flying) +
    (1 | bat_roosting) + (1 | group), data = d, ziformula = ~1, family = nbinom2)
check_overdispersion(t)
# Overdispersion test

       dispersion ratio =   0.582
  Pearson's Chi-Squared = 255.097
                p-value =       1
Code
check_zeroinflation(t)
# Check for zero-inflation

   Observed zeros: 191
  Predicted zeros: 30
            Ratio: 0.16
Code
t2 <- tidy(t, conf.int = TRUE, exponentiate = T, effects = "fixed",
    conf.method = "profile")
t2
# A tibble: 4 × 9
  effect component term  estimate std.error statistic p.value conf.low conf.high
  <chr>  <chr>     <chr>    <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 fixed  cond      (Int…    3.00     1.60        2.06 3.95e-2     1.08      8.91
2 fixed  cond      flig…    1.58     0.173       4.18 2.98e-5     1.26      1.95
3 fixed  cond      log_…    1.45     0.188       2.83 4.62e-3     1.12      1.87
4 fixed  zi        (Int…    0.644    0.0722     -3.93 8.65e-5     1.08      8.91
Code
# for every minute of flight time, the response call count
# increases by a factor of 1.58 [1.26, 1.95]

##### Effect of inquiry calling on response calling------------

# plot with log counts
d %>%
    ggplot(aes(x = log_inquiry, y = log(n_response + 1))) + geom_point(size = 2) +
    geom_smooth(method = "lm") + xlab("log inquiry call count") +
    ylab("log (x+1) response call count")

Code
# plot negative binomial curve plot all bats
(t1 <- d %>%
    mutate(label = "all pairs") %>%
    mutate(escape = as.logical(escape)) %>%
    ggplot(aes(x = log_inquiry, y = n_response)) + facet_wrap(~label) +
    geom_point(size = 2, alpha = 0.5) + geom_smooth(method = "glm.nb") +
    xlab("log of inquiry call count") + ylab("response call count")) +
    theme(legend.position = "none")

Code
# plot with kinship
(t2 <- d %>%
    filter(!is.na(kinship)) %>%
    mutate(escape = as.logical(escape)) %>%
    mutate(kin = if_else(kinship > 0.1, "kin", "nonkin")) %>%
    mutate(kin = factor(kin, levels = c("nonkin", "kin"))) %>%
    ggplot(aes(x = log_inquiry, y = (n_response), color = kin, group = kin)) +
    facet_wrap(~kin) + geom_point(size = 2, alpha = 0.5, aes(shape = escape)) +
    geom_smooth(method = "glm.nb") + xlab("log of inquiry call count") +
    ylab("response call count") + scale_color_manual(values = c("darkgrey",
    "darkred")) + theme(legend.position = "none") + theme(axis.text.y = element_blank(),
    axis.title.y = element_blank(), legend.position = "none"))

Code
# plot as linear model
(t3 <- d %>%
    filter(!is.na(kinship)) %>%
    mutate(kin = if_else(kinship > 0.1, "kin", "nonkin")) %>%
    mutate(kin = factor(kin, levels = c("nonkin", "kin"))) %>%
    ggplot(aes(x = log_inquiry, y = log(n_response + 1), color = kin,
        group = kin)) + facet_wrap(~kin) + geom_point(size = 2, alpha = 0.5) +
    geom_smooth(method = "lm") + xlab("log of inquiry call count") +
    ylab("log of response call count") + scale_color_manual(values = c("darkgrey",
    "darkred")) + theme(legend.position = "none"))

Code
(inquiry.response <- t1 + t2 + plot_layout(ncol = 2) + plot_layout(widths = c(1,
    2)))

Code
# save as PDF ggsave( 'inquiry_response.pdf', plot =
# inquiry.response, width = 14, height = 7, units = 'in', dpi =
# 300)

# for plots, get residuals from negative binomial mixed effect
# model predicting number of responses by number of inquiry
# calls
fit <- glmmTMB(n_response ~ log_inquiry + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~1, family = nbinom2)
# get difference between observed response count and predicted
# response count add this response score to data for plotting
# this value represents response calling more than expected from
# inquiry calls
d$resid_response <- d$n_response - predict(fit, type = "response",
    newdata = d, allow.new.levels = T)


# plot number of inquiry calls by kinship with partner
d %>%
    filter(!is.na(kinship)) %>%
    mutate(kin = as.character(kinship)) %>%
    ggplot(aes(x = kin, y = n_inquiry, group = kin)) + geom_violin(fill = NA,
    width = 1) + geom_jitter(size = 2, width = 0.1, height = 0, aes(color = udyad.sex)) +
    xlab("kinship") + ylab("number of inquiry calls") + theme(legend.position = "top",
    legend.title = element_blank())

Code
# plot number of response calls by kinship with partner
d %>%
    filter(!is.na(kinship)) %>%
    mutate(kin = as.character(kinship)) %>%
    ggplot(aes(x = kin, y = n_response, group = kin)) + geom_violin(fill = NA,
    width = 1) + geom_jitter(size = 2, width = 0.1, height = 0, aes(color = udyad.sex)) +
    xlab("kinship") + ylab("number of response calls") + theme(legend.position = "top",
    legend.title = element_blank())

Code
# plot residual response call variation
d %>%
    filter(!is.na(kinship)) %>%
    mutate(kin = as.character(kinship)) %>%
    ggplot(aes(x = kin, y = resid_response, group = kin)) + geom_violin(fill = NA,
    width = 1) + geom_jitter(size = 2, width = 0.1, height = 0, aes(color = udyad.sex)) +
    xlab("kinship") + ylab("residual response call variation") + theme(legend.position = "top",
    legend.title = element_blank())

Code
# create function to plot variable by kinship with means and 95%
# CIs
plot_kinship_effect <- function(d = d, y = d$sri, label = "label") {
    # plot association by kinship
    set.seed(123)
    means <- d %>%
        mutate(y = y) %>%
        filter(!is.na(kinship)) %>%
        group_by(dyad) %>%
        summarize(kinship = mean(kinship), y = mean(y), udyad.sex = first(udyad.sex)) %>%
        mutate(kin = ifelse(kinship > 0, "kin", "nonkin")) %>%
        dplyr::select(kin, y) %>%
        boot_ci2(y = .$y, x = .$kin, bca = T) %>%
        rename(kin = effect) %>%
        mutate(kin = factor(kin, levels = c("nonkin", "kin")))
    points <- d %>%
        mutate(y = y) %>%
        filter(!is.na(kinship)) %>%
        group_by(dyad) %>%
        summarize(kinship = mean(kinship), y = mean(y), udyad.sex = first(udyad.sex)) %>%
        mutate(kin = ifelse(kinship > 0, "kin", "nonkin")) %>%
        mutate(kin = factor(kin, levels = c("nonkin", "kin")))
    # plot means and 95% CI
    means %>%
        ggplot(aes(x = kin, y = mean, color = kin)) + geom_jitter(data = points,
        aes(y = y), size = 1, alpha = 0.5, width = 0.1) + geom_point(position = position_nudge(x = 0.25),
        size = 3) + geom_errorbar(aes(ymin = low, ymax = high, width = 0.1),
        position = position_nudge(x = 0.25), size = 1) + xlab("") +
        ylab(label) + scale_color_manual(values = c("darkgrey", "darkred")) +
        coord_flip() + theme_bw() + theme(legend.position = "none")
}

# plot effect of kinship on SRI
(p1 <- plot_kinship_effect(d, d$sri, label = "association (SRI)"))

Code
# plot effect of kinship on inquiry calling
(p2 <- plot_kinship_effect(d, d$n_inquiry/(d$flight_time/60), label = "inquiry calls per min of flight time"))

Code
# plot effect of kinship on response calling
(p3 <- plot_kinship_effect(d, d$n_response/(d$flight_time/60), label = "response calls per min of flight time"))

Code
# plot effect of kinship on response calling
(p4 <- plot_kinship_effect(d, d$resid_response, label = "residual response call variation"))

Code
# plot all
(kinship.effects <- p1 + p2 + p3 + p4 + plot_layout(ncol = 1))

Code
# save as PDF ggsave( 'kinship_effects.pdf', plot =
# kinship.effects, scale = 1, width = 3.5, height = 7, units =
# 'in', dpi = 300)
Code
# plot number of inquiry calls by association and kinship (by
# dyad type)
t1 <- d %>%
    filter(!is.na(kinship)) %>%
    mutate(kin = ifelse(kinship > 0.1, "kin", "nonkin")) %>%
    ggplot(aes(x = sri, y = log_inquiry)) + facet_wrap(~dyad.sex,
    scales = "free_y") + geom_point(aes(color = kin), size = 2, alpha = 0.6) +
    geom_smooth(method = "lm") + xlab("association rate (simple ratio index)") +
    ylab("inquiry call count (log transformed)") + scale_colour_manual(values = c("darkred",
    "grey")) + theme_bw() + theme(legend.position = "none", legend.title = element_blank())

# plot number of response calls by association and kinship (by
# dyad type)
t2 <- d %>%
    filter(!is.na(kinship)) %>%
    mutate(kin = ifelse(kinship > 0.1, "kin", "nonkin")) %>%
    ggplot(aes(x = sri, y = log(n_response + 1))) + facet_wrap(~dyad.sex,
    scales = "free_y") + geom_point(aes(color = kin), size = 2, alpha = 0.6) +
    geom_smooth(method = "lm") + xlab("association rate (simple ratio index)") +
    ylab("response call count (log transformed)") + scale_colour_manual(values = c("darkred",
    "grey")) + theme_bw() + theme(legend.position = "none", legend.title = element_blank())

# plot residual response calling by association (by dyad type)-
# controlling for inquiry calls
t3 <- d %>%
    filter(!is.na(kinship)) %>%
    mutate(kin = ifelse(kinship > 0.1, "kin", "nonkin")) %>%
    ggplot(aes(x = sri, y = resid_response)) + facet_wrap(~dyad.sex,
    scales = "free_y") + geom_point(aes(color = kin), size = 2, alpha = 0.6) +
    geom_smooth(method = "lm") + xlab("association rate (simple ratio index)") +
    ylab("residual response call variation") + scale_colour_manual(values = c("darkred",
    "grey")) + theme_bw() + theme(legend.position = "none", legend.title = element_blank())

# combine plots
(by.sex.plot <- t1 + t2 + t3 + plot_layout(ncol = 1) + plot_annotation(tag_levels = "A"))

Code
# save as PDF ggsave( 'by_sex_plot.pdf', plot = by.sex.plot,
# scale = 1, width = 7, height = 14, units = 'in', dpi = 300)

3.3.1 Effect of kinship and association on response calling

Code
# poisson model is overdispersed
t <- glmer(n_response ~ kinship * sri + log_inquiry + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), family = poisson,
    data = d)
check_overdispersion(t)
# Overdispersion test

       dispersion ratio =    72.342
  Pearson's Chi-Squared = 22570.653
                p-value =   < 0.001
Code
# fit negative binomial model (NBM)
t <- glmmTMB(n_response ~ kinship * sri + log_inquiry + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~0, family = nbinom2)
check_overdispersion(t)
# Overdispersion test

       dispersion ratio =   0.335
  Pearson's Chi-Squared = 104.319
                p-value =       1
Code
check_zeroinflation(t, tolerance = 0.05)
# Check for zero-inflation

   Observed zeros: 133
  Predicted zeros: 116
            Ratio: 0.87
Code
AIC(t)  # 2700
[1] 2699.808
Code
BIC(t)  # 2734
[1] 2733.723
Code
# fit zero-inflated NBM
fit <- glmmTMB(n_response ~ kinship * sri + log_inquiry + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~1, family = nbinom2)
AIC(fit)  #2649
[1] 2648.867
Code
BIC(fit)  # 2687
[1] 2686.55
Code
summary(fit)
 Family: nbinom2  ( log )
Formula:          
n_response ~ kinship * sri + log_inquiry + offset(log(flight_time)) +  
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation:              ~1
Data: d

     AIC      BIC   logLik deviance df.resid 
  2648.9   2686.6  -1314.4   2628.9      310 

Random effects:

Conditional model:
 Groups       Name        Variance  Std.Dev. 
 bat_flying   (Intercept) 3.792e-08 0.0001947
 bat_roosting (Intercept) 2.296e-01 0.4791231
 group        (Intercept) 1.678e-01 0.4096551
Number of obs: 320, groups:  bat_flying, 50; bat_roosting, 50; group, 21

Dispersion parameter for nbinom2 family (): 0.707 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -1.2182     0.6294  -1.935   0.0529 .
kinship      -3.3360     1.9159  -1.741   0.0817 .
sri          -0.7478     0.6424  -1.164   0.2443  
log_inquiry   0.1897     0.1154   1.644   0.1002  
kinship:sri   5.2163     2.7487   1.898   0.0577 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.4362     0.1242  -3.512 0.000444 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# fit zero-inflated NBM without interaction
fit <- glmmTMB(n_response ~ kinship + sri + log_inquiry + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~1, family = nbinom2)
AIC(fit)  #2650
[1] 2650.47
Code
BIC(fit)  # 2684
[1] 2684.385
Code
summary(fit)
 Family: nbinom2  ( log )
Formula:          
n_response ~ kinship + sri + log_inquiry + offset(log(flight_time)) +  
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation:              ~1
Data: d

     AIC      BIC   logLik deviance df.resid 
  2650.5   2684.4  -1316.2   2632.5      311 

Random effects:

Conditional model:
 Groups       Name        Variance  Std.Dev. 
 bat_flying   (Intercept) 4.139e-08 0.0002035
 bat_roosting (Intercept) 1.976e-01 0.4445021
 group        (Intercept) 1.687e-01 0.4107856
Number of obs: 320, groups:  bat_flying, 50; bat_roosting, 50; group, 21

Dispersion parameter for nbinom2 family (): 0.678 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -1.6725     0.5816  -2.876  0.00403 **
kinship       0.2071     0.5173   0.400  0.68894   
sri          -0.2454     0.5926  -0.414  0.67882   
log_inquiry   0.2338     0.1138   2.054  0.03995 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.4456     0.1256  -3.548 0.000388 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# get predicted values
d$n_response_predicted <- predict(fit, type = "response", newdata = d,
    allow.new.levels = T)

# plot model performance
d %>%
    filter(!is.na(kinship)) %>%
    mutate(kinship = kinship > 0) %>%
    ggplot(aes(x = n_response, y = n_response_predicted)) + geom_point(size = 2,
    aes(color = kinship)) + geom_smooth(method = "lm") + xlab("observed count of response calls") +
    ylab("predicted count of response calls") + scale_color_manual(values = c("darkgrey",
    "darkred")) + theme(legend.position = "top")

Code
# get relationship between predicted and observed
summary(lm(n_response_predicted ~ n_response, data = d))  # r-squared= 0.46

Call:
lm(formula = n_response_predicted ~ n_response, data = d)

Residuals:
     Min       1Q   Median       3Q      Max 
-105.210  -18.559   -1.269   18.411  109.081 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 58.98540    2.02304   29.16   <2e-16 ***
n_response   0.18176    0.01096   16.58   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 31.09 on 318 degrees of freedom
  (126 observations deleted due to missingness)
Multiple R-squared:  0.4638,    Adjusted R-squared:  0.4621 
F-statistic:   275 on 1 and 318 DF,  p-value: < 2.2e-16
Code
# get model coefficients with 95% CIs
coefs <- tidy(fit, conf.int = TRUE, exponentiate = F, effects = "fixed",
    conf.method = "profile") %>%
    filter(component == "cond") %>%
    dplyr::select(-effect, -component) %>%
    mutate(type = "full model")
# coefs

# save model results write.csv(coefs,
# file='response_model_results.csv')

# effect of SRI with kin only
fitk <- glmmTMB(n_response ~ sri + log_inquiry + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d %>%
    filter(kinship > 0), ziformula = ~1, family = nbinom2)
coefs.k <- tidy(fitk, conf.int = TRUE, exponentiate = F, effects = "fixed",
    conf.method = "profile") %>%
    filter(component == "cond") %>%
    dplyr::select(-effect, -component) %>%
    mutate(type = "kin only")
# coefs.k

# effect of SRI with nonkin only
fitn <- glmmTMB(n_response ~ sri + log_inquiry + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d %>%
    filter(kinship == 0), ziformula = ~1, family = nbinom2)
# get 95% CI
coefs.n <- tidy(fitn, conf.int = TRUE, exponentiate = F, effects = "fixed",
    conf.method = "profile") %>%
    filter(component == "cond") %>%
    dplyr::select(-effect, -component) %>%
    mutate(type = "nonkin only")
# coefs.n

# plot model results
theme_set(theme_bw(base_size = 12))
plot1 <- coefs %>%
    filter(term != "(Intercept)") %>%
    mutate(term = case_when(term == "sri" ~ "association", term ==
        "kinship" ~ "kinship", term == "log_inquiry" ~ "log count of inquiry calls",
        term == "kinship:sri" ~ "association x kinship interaction")) %>%
    mutate(term = factor(term, levels = c("log count of inquiry calls",
        "association x kinship interaction", "association", "kinship"))) %>%
    mutate(term = fct_rev(term)) %>%
    ggplot(aes(x = estimate, y = term)) + geom_point(size = 2) + geom_errorbarh(aes(xmin = conf.low,
    xmax = conf.high, height = 0.2), size = 1) + geom_vline(xintercept = 0,
    linetype = "dashed") + ylab("") + xlab("coefficient") + coord_cartesian(xlim = c(-1.5,
    1.5)) + theme(axis.text = element_text(size = 12), strip.text = element_text(size = 12))
plot1

Code
# save as PDF ggsave( 'response_model_results.pdf', plot =
# plot1, scale = 1, width = 5, height = 2, units = 'in', dpi =
# 300)


# check that sri and kinship do not predict response calls as
# single predictors fit zero-inflated negative binomial model
# with only sri
t <- glmmTMB(n_response ~ sri + log_inquiry + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~1, family = nbinom2)
summary(t)
 Family: nbinom2  ( log )
Formula:          
n_response ~ sri + log_inquiry + offset(log(flight_time)) + (1 |  
    bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation:              ~1
Data: d

     AIC      BIC   logLik deviance df.resid 
  3579.0   3611.6  -1781.5   3563.0      428 

Random effects:

Conditional model:
 Groups       Name        Variance  Std.Dev. 
 bat_flying   (Intercept) 2.229e-08 0.0001493
 bat_roosting (Intercept) 6.142e-01 0.7836871
 group        (Intercept) 5.587e-02 0.2363637
Number of obs: 436, groups:  bat_flying, 72; bat_roosting, 73; group, 23

Dispersion parameter for nbinom2 family (): 0.646 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.1155     0.5679  -3.725 0.000195 ***
sri          -0.1597     0.4814  -0.332 0.740114    
log_inquiry   0.3257     0.1102   2.955 0.003122 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.4676     0.1155  -4.048 5.17e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# fit zero-inflated negative binomial model with only kinship
t <- glmmTMB(n_response ~ kinship + log_inquiry + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~1, family = nbinom2)
summary(t)
 Family: nbinom2  ( log )
Formula:          
n_response ~ kinship + log_inquiry + offset(log(flight_time)) +  
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation:              ~1
Data: d

     AIC      BIC   logLik deviance df.resid 
  2648.6   2678.8  -1316.3   2632.6      312 

Random effects:

Conditional model:
 Groups       Name        Variance  Std.Dev. 
 bat_flying   (Intercept) 2.750e-08 0.0001658
 bat_roosting (Intercept) 2.023e-01 0.4498245
 group        (Intercept) 1.786e-01 0.4226129
Number of obs: 320, groups:  bat_flying, 50; bat_roosting, 50; group, 21

Dispersion parameter for nbinom2 family (): 0.681 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.7899     0.5084  -3.521  0.00043 ***
kinship       0.1260     0.4792   0.263  0.79257    
log_inquiry   0.2263     0.1125   2.011  0.04429 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.4447     0.1254  -3.545 0.000393 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
######## Effect of kinship and association on inquiry
######## calling---------

# poisson model is overdispersed
t <- glmer(n_inquiry ~ kinship * sri + log(n_response + 1) + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), family = poisson,
    data = d)
check_overdispersion(t)
# Overdispersion test

       dispersion ratio =    6.731
  Pearson's Chi-Squared = 2100.164
                p-value =  < 0.001
Code
# fit negative binomial model (NBM) with log responses
t <- glmmTMB(n_inquiry ~ kinship * sri + log(n_response + 1) + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~0, family = nbinom2)
check_overdispersion(t)
# Overdispersion test

       dispersion ratio =   0.820
  Pearson's Chi-Squared = 255.132
                p-value =   0.991
Code
AIC(t)  # 3130
[1] 3129.521
Code
BIC(t)  # 3163
[1] 3163.436
Code
# fit NBM with responses
t <- glmmTMB(n_inquiry ~ kinship * sri + n_response + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~0, family = nbinom2)
check_overdispersion(t)
# Overdispersion test

       dispersion ratio =   0.813
  Pearson's Chi-Squared = 252.950
                p-value =   0.993
Code
AIC(t)  #3130
[1] 3130.339
Code
BIC(t)  # 3164
[1] 3164.254
Code
summary(t)
 Family: nbinom2  ( log )
Formula:          
n_inquiry ~ kinship * sri + n_response + offset(log(flight_time)) +  
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d

     AIC      BIC   logLik deviance df.resid 
  3130.3   3164.3  -1556.2   3112.3      311 

Random effects:

Conditional model:
 Groups       Name        Variance  Std.Dev. 
 bat_flying   (Intercept) 1.833e-01 4.282e-01
 bat_roosting (Intercept) 2.748e-09 5.242e-05
 group        (Intercept) 3.400e-02 1.844e-01
Number of obs: 320, groups:  bat_flying, 50; bat_roosting, 50; group, 21

Dispersion parameter for nbinom2 family (): 4.68 

Conditional model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.3629986  0.1517410  -8.982   <2e-16 ***
kinship     -0.5338331  0.6126168  -0.871    0.384    
sri         -0.1397377  0.1970890  -0.709    0.478    
n_response   0.0001023  0.0002075   0.493    0.622    
kinship:sri  0.4310938  0.9048669   0.476    0.634    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# fit NBM with responses without interaction
fit2 <- glmmTMB(n_inquiry ~ kinship + sri + n_response + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~0, family = nbinom2)
AIC(fit2)  #3129
[1] 3128.567
Code
BIC(fit2)  # 3159
[1] 3158.713
Code
summary(fit2)
 Family: nbinom2  ( log )
Formula:          
n_inquiry ~ kinship + sri + n_response + offset(log(flight_time)) +  
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d

     AIC      BIC   logLik deviance df.resid 
  3128.6   3158.7  -1556.3   3112.6      312 

Random effects:

Conditional model:
 Groups       Name        Variance  Std.Dev. 
 bat_flying   (Intercept) 1.861e-01 0.4314410
 bat_roosting (Intercept) 2.571e-09 0.0000507
 group        (Intercept) 3.441e-02 0.1855110
Number of obs: 320, groups:  bat_flying, 50; bat_roosting, 50; group, 21

Dispersion parameter for nbinom2 family (): 4.69 

Conditional model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.3851696  0.1446519  -9.576   <2e-16 ***
kinship     -0.2524912  0.1637266  -1.542    0.123    
sri         -0.1002280  0.1789141  -0.560    0.575    
n_response   0.0001123  0.0002065   0.544    0.587    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# get predicted values
d$n_inquiry_predicted <- predict(fit2, type = "response", newdata = d,
    allow.new.levels = T)

# plot model performance
d %>%
    filter(!is.na(kinship)) %>%
    mutate(kinship = kinship > 0) %>%
    ggplot(aes(x = n_inquiry, y = n_inquiry_predicted)) + geom_point(size = 2,
    aes(color = kinship)) + geom_smooth(method = "lm") + xlab("observed count of inquiry calls") +
    ylab("predicted count of inquiry calls") + scale_color_manual(values = c("darkgrey",
    "darkred")) + theme(legend.position = "top")

Code
# get relationship between predicted and observed
summary(lm(n_inquiry_predicted ~ n_inquiry, data = d))  # r-squared= 0.71

Call:
lm(formula = n_inquiry_predicted ~ n_inquiry, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-92.258 -12.967  -0.546  11.381  87.999 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 19.40700    2.06189   9.412   <2e-16 ***
n_inquiry    0.71686    0.02554  28.063   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 19.59 on 318 degrees of freedom
  (126 observations deleted due to missingness)
Multiple R-squared:  0.7124,    Adjusted R-squared:  0.7114 
F-statistic: 787.5 on 1 and 318 DF,  p-value: < 2.2e-16
Code
# get model coefficients
summary(fit2)
 Family: nbinom2  ( log )
Formula:          
n_inquiry ~ kinship + sri + n_response + offset(log(flight_time)) +  
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d

     AIC      BIC   logLik deviance df.resid 
  3128.6   3158.7  -1556.3   3112.6      312 

Random effects:

Conditional model:
 Groups       Name        Variance  Std.Dev. 
 bat_flying   (Intercept) 1.861e-01 0.4314410
 bat_roosting (Intercept) 2.571e-09 0.0000507
 group        (Intercept) 3.441e-02 0.1855110
Number of obs: 320, groups:  bat_flying, 50; bat_roosting, 50; group, 21

Dispersion parameter for nbinom2 family (): 4.69 

Conditional model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.3851696  0.1446519  -9.576   <2e-16 ***
kinship     -0.2524912  0.1637266  -1.542    0.123    
sri         -0.1002280  0.1789141  -0.560    0.575    
n_response   0.0001123  0.0002065   0.544    0.587    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
coefs2 <- tidy(fit2, conf.int = TRUE, exponentiate = F, effects = "fixed",
    conf.method = "profile") %>%
    filter(component == "cond") %>%
    dplyr::select(-effect, -component) %>%
    mutate(type = "full model")
# coefs2

# save model results write.csv(coefs2,
# file='inquiry_model_results.csv')

# plot model results
theme_set(theme_bw(base_size = 12))
plot2 <- coefs2 %>%
    filter(term != "(Intercept)") %>%
    mutate(term = case_when(term == "sri" ~ "association", term ==
        "kinship" ~ "kinship", term == "n_response" ~ "response calls")) %>%
    mutate(term = factor(term, levels = c("response calls", "association",
        "kinship"))) %>%
    mutate(term = fct_rev(term)) %>%
    ggplot(aes(x = estimate, y = term)) + geom_point(size = 2) + geom_errorbarh(aes(xmin = conf.low,
    xmax = conf.high, height = 0.2), size = 1) + geom_vline(xintercept = 0,
    linetype = "dashed") + ylab("") + xlab("estimate (log odds)") +
    theme(axis.text = element_text(size = 12), strip.text = element_text(size = 12))
plot2

Code
# save as PDF ggsave( 'inquiry_model_results.pdf', plot = plot1,
# scale = 1, width = 7, height = 3.5, units = 'in', dpi = 300)

# get relative amount of variance explained by random intercept
tidy(fit2, conf.int = F, exponentiate = F) %>%
    filter(effect == "ran_pars") %>%
    mutate(variance = estimate^2) %>%
    dplyr::select(effect, group, variance) %>%
    group_by(effect) %>%
    mutate(sum = sum(variance)) %>%
    mutate(prop = variance/sum)
# A tibble: 3 × 5
# Groups:   effect [1]
  effect   group             variance   sum         prop
  <chr>    <chr>                <dbl> <dbl>        <dbl>
1 ran_pars bat_flying   0.186         0.221 0.844       
2 ran_pars bat_roosting 0.00000000257 0.221 0.0000000117
3 ran_pars group        0.0344        0.221 0.156       
Code
# confirm no effect with single predictors association
t <- glmmTMB(n_inquiry ~ sri + n_response + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~0, family = nbinom2)
summary(t)
 Family: nbinom2  ( log )
Formula:          
n_inquiry ~ sri + n_response + offset(log(flight_time)) + (1 |  
    bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d

     AIC      BIC   logLik deviance df.resid 
  4257.0   4285.5  -2121.5   4243.0      429 

Random effects:

Conditional model:
 Groups       Name        Variance Std.Dev.
 bat_flying   (Intercept) 0.23098  0.48060 
 bat_roosting (Intercept) 0.00448  0.06694 
 group        (Intercept) 0.01246  0.11164 
Number of obs: 436, groups:  bat_flying, 72; bat_roosting, 73; group, 23

Dispersion parameter for nbinom2 family ():  5.1 

Conditional model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.3782473  0.1046922 -13.165   <2e-16 ***
sri         -0.1730846  0.1317072  -1.314   0.1888    
n_response   0.0003027  0.0001691   1.790   0.0734 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# kinship
t <- glmmTMB(n_inquiry ~ kinship + n_response + offset(log(flight_time)) +
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
    ziformula = ~0, family = nbinom2)
summary(t)
 Family: nbinom2  ( log )
Formula:          
n_inquiry ~ kinship + n_response + offset(log(flight_time)) +  
    (1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d

     AIC      BIC   logLik deviance df.resid 
  3126.9   3153.3  -1556.4   3112.9      313 

Random effects:

Conditional model:
 Groups       Name        Variance  Std.Dev. 
 bat_flying   (Intercept) 1.825e-01 4.272e-01
 bat_roosting (Intercept) 3.244e-09 5.695e-05
 group        (Intercept) 3.491e-02 1.868e-01
Number of obs: 320, groups:  bat_flying, 50; bat_roosting, 50; group, 21

Dispersion parameter for nbinom2 family (): 4.67 

Conditional model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.4480670  0.0905664 -15.989   <2e-16 ***
kinship     -0.2738852  0.1594311  -1.718   0.0858 .  
n_response   0.0001283  0.0002052   0.625   0.5319    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Takeaways

  • Individual ID but not kinship or association explain calling rates

 


 

Session information

R version 4.2.2 Patched (2022-11-10 r83330)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=es_ES.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=es_ES.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=es_ES.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] MASS_7.3-58.2         performance_0.10.3    patchwork_1.1.2      
 [4] boot_1.3-28           broom.mixed_0.2.9.4   glmmTMB_1.1.7        
 [7] lme4_1.1-33           Matrix_1.5-1          mapproj_1.2.11       
[10] maps_3.4.1            geosphere_1.5-18      lubridate_1.9.2      
[13] forcats_1.0.0         stringr_1.5.0         dplyr_1.1.0          
[16] purrr_1.0.1           readr_2.1.4           tidyr_1.3.0          
[19] tibble_3.2.1          tidyverse_2.0.0       RColorBrewer_1.1-3   
[22] statnet_2019.6        tsna_0.3.5            ergm.count_4.1.1     
[25] tergm_4.2.0           networkDynamic_0.11.3 ergm_4.5.0           
[28] igraphdata_1.0.1      ecodist_2.0.9         ggmap_3.0.1          
[31] ggplot2_3.4.3         assortnet_0.20        asnipe_1.1.17        
[34] bipartite_2.18        sna_2.7-1             network_1.18.1       
[37] statnet.common_4.9.0  vegan_2.6-4           lattice_0.20-45      
[40] permute_0.9-7         igraph_1.5.1         

loaded via a namespace (and not attached):
 [1] minqa_1.2.5              colorspace_2.1-0         ellipsis_0.3.2          
 [4] rprojroot_2.0.3          estimability_1.4.1       xaringanExtra_0.7.0     
 [7] rstudioapi_0.14          farver_2.1.1             listenv_0.8.0           
[10] furrr_0.3.1              remotes_2.4.2            mvtnorm_1.1-3           
[13] fansi_1.0.4              codetools_0.2-19         splines_4.2.2           
[16] cachem_1.0.8             robustbase_0.99-0        knitr_1.43              
[19] spam_2.9-1               jsonlite_1.8.7           nloptr_2.0.3            
[22] packrat_0.8.1            broom_1.0.4              cluster_2.1.4           
[25] png_0.1-8                compiler_4.2.2           httr_1.4.6              
[28] backports_1.4.1          emmeans_1.8.6            fastmap_1.1.1           
[31] cli_3.6.1                formatR_1.12             htmltools_0.5.6         
[34] tools_4.2.2              dotCall64_1.0-2          coda_0.19-4             
[37] gtable_0.3.4             glue_1.6.2               Rcpp_1.0.11             
[40] rle_0.9.2                vctrs_0.6.3              nlme_3.1-162            
[43] insight_0.19.2           xfun_0.40                globals_0.16.1          
[46] rbibutils_2.2.15         trust_0.1-8              timechange_0.2.0        
[49] lifecycle_1.0.3          future_1.28.0            DEoptimR_1.1-2          
[52] scales_1.2.1             hms_1.1.2                parallel_4.2.2          
[55] ergm.multi_0.2.0         TMB_1.9.6                fields_15.2             
[58] yaml_2.3.7               memoise_2.0.1            stringi_1.7.12          
[61] RgoogleMaps_1.4.5.3      Rdpack_2.5               rlang_1.1.1             
[64] pkgconfig_2.0.3          bitops_1.0-7             evaluate_0.21           
[67] lpSolveAPI_5.5.2.0-17.10 labeling_0.4.3           htmlwidgets_1.5.4       
[70] tidyselect_1.2.0         parallelly_1.32.1        plyr_1.8.7              
[73] magrittr_2.0.3           R6_2.5.1                 generics_0.1.3          
[76] pillar_1.9.0             withr_2.5.0              mgcv_1.8-41             
[79] sp_2.0-0                 crayon_1.5.2             utf8_1.2.3              
[82] networkLite_1.0.5        tzdb_0.3.0               rmarkdown_2.21          
[85] jpeg_0.1-9               grid_4.2.2               sketchy_1.0.2           
[88] digest_0.6.33            xtable_1.8-4             numDeriv_2016.8-1.1     
[91] munsell_0.5.0            viridisLite_0.4.2